home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 22
/
Cream of the Crop 22.iso
/
program
/
tn2.zip
/
GAUSSIAN.T
< prev
next >
Wrap
Text File
|
1996-11-15
|
3KB
|
163 lines
%
% "gaussian.t" performs gaussian elimination on
% an N x N matrix and solves for x in:
%
% A x = b
%
% Sample program for the T Interpreter by:
%
% Stephen R. Schmitt
% 962 Depot Road
% Boxborough, MA 01719
%
const N : int := 3
var a : array[N,N+1] of real
var x : array[N] of real
program
var i, j : int
a[0,0] := 1 % a00
a[0,1] := 3 % a01
a[0,2] := -4 % a02
a[0,3] := 8 % b0
a[1,0] := 1 % a10
a[1,1] := 1 % a11
a[1,2] := -2 % a12
a[1,3] := 2 % b1
a[2,0] := -1 % a20
a[2,1] := -2 % a21
a[2,2] := 5 % a22
a[2,3] := -1 % b2
x[0] := 0 % initialize x
x[1] := 0
x[2] := 0
if eliminate then % lower diagonal elements
substitute % and solve for x
put "solution:"
put ""
for i := 0 ... N-1 do
put " x[", i, "] = ", x[i]
end for
else
put "singular matrix!"
end if
end program
%
% "eliminate" converts A to upper diagonal form
%
% returns: nothing
%
function eliminate : boolean
var i, j, k, max_row : int
var t : real
var ok : boolean := true
for i := 0 ... N-1 do
%
% find row with maximum in column i
%
max_row := i
for j := i+1 ... N-1 do
if abs( a[j,i] ) > abs( a[max_row,i] ) then
max_row := j
end if
end for
%
% swap max row with row i of A and b
%
for k := i ... N do
t := a[i,k]
a[i,k] := a[max_row,k]
a[max_row,k] := t
end for
%
% eliminate lower diagonal elements
%
for j := i+1 ... N-1 do
for decreasing k := N ... i do
if a[i,i] = 0.0 then
ok := false
else
a[j,k] := a[j,k] - a[i,k] * a[j,i] / a[i,i]
end if
end for
end for
end for
return ok
end function
%
% "substitute" computes the values of x starting from the bottom
%
% returns: nothing
%
procedure substitute
var j, k : int
var t : real
for decreasing j := N-1 ... 0 do
t := 0.0
for k := j+1 ... N-1 do
t := t + a[j,k] * x[k]
end for
x[j] := ( a[j,N] - t ) / a[j,j]
end for
end procedure
%
% "abs" returns the absolute value of a real argument
%
% returns: real number
%
function abs( x : real ) : real
if x < 0 then
x := -x
end if
return x
end function